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ABSTRACT 

Context. In massive binaries, the powerful stellar winds of the two stars collide, leading to the formation of shock-dominated envi- 
ronments that can be modelled only in 3D. 

Aims. We investigate the morphology of the collision front between the stellar winds of binary components in two long-period binary 
systems, one consisting of a hydrogen rich Wolf-Rayet star (WNL) and an O-star and the other of a Luminous Blue Variable (LBV) 
and an O-star. Specifically, we follow the development and evolution of instabilities that form in such a shell, if it is sufficiently 
compressed, due to both the wind interaction and the orbital motion. 

Methods. We use MPI-AMRVAC to time-integrate the equations of hydrodynamics, combined with optically thin radiative cooling, 
on an adaptive mesh 3D grid. Using parameters for generic binary systems, we simulate the interaction between the winds of the two 
stars. 

Results. The WNL + O star binary shows a typical example of an adiabatic wind collision. The resulting shell is thick and smooth, 
showing no instabilities. On the other hand, the shell created by the collision of the O star wind with the LBV wind, combined with 
the orbital motion of the binary components, is susceptible to thin shell instabilities, which create a highly structured morphology. 
We identify the nature of the instabilities as both linear and non-linear thin-shell instabilities, with distinct differences between the 
leading and the trailing parts of the collision front. We also find that for binaries containing a star with a (relatively) slow wind, the 
global shape of the shell is determined more by the slow wind velocity and the orbital motion of the binary, than the ram pressure 
balance between the two winds. 

Conclusions. The interaction between massive binary winds needs further parametric exploration, to identify the role and dynamical 
importance of multiple instabilities at the collision front, as shown here for an LBV + O star system. 
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1. Modelling circumstellar environments 

Massive stars lose a significant amount of their mass during their 
evolution. The mass is lost in the form of both quasi-continuous 
winds and sporadic eruptions. Changes in mass-loss rate and ve- 
locity create a series of hydrodynamical interactions, shaping the 
circumstellar medium. The result of these interactions can be ob- 
served in the form of circumstellar nebulae, which are temporary 
high density structures surrounding the star. 

The formation of such structures around single stars 
has been successfully modelled in quite some de tai l (e. 



Frgver et all I2003L l200i 
However, since most 



Garcia -Seeur a et alJ TT996bllal ; 
van Marie et alJ I2005L I2007L l2Qq 
massive stars occur in binaries, the value of such single star 
simulations is limited. At the same time, single star circumstellar 
bubbles can be studied in restricted dimensionality, allowing 
very high resolution parametric studies, highlighting the role of 
various thin- shell instabilities in circumstellar morphology. For 
a binary, even the symmetric case of two identical stars requires 
3D computations, at considerably increased computational 
costs. 

Recent work exploitin g Smooth Particle Hydrodynamics 
(SPH) methodologies (e.g. lOkazaki et al.ll2008l) investigated the 
3D wind- wind interactions believed to occur in the Eta Carinae 
star system. However, standard SPH algorithms face difficulties 



in the vicinity of strong shocks, while the collision between two 
strong stellar winds inherently leads to their creation. Modelling 
the wind collision using grid-based methods has usually fo- 
cused on the details of the inter action front in between the two 
stars, such as in the work by IWaldei] (Il998l) ; iFolini & Wa lder 
( 2000). More recently, models emerged where both symmetric 
and asymmetric systems of massi ve O-stars, on circular or ellip- 
tic orb i ts, were follow ed in detail (|Pittard 2009l; |Pittard & Parkinl 
|2010|) . iPittardl d2009h in particular uses Riemann- solver based 
discretizations, on a fixed 3D grid, solving hydrodynamics to- 
gether with optically thin radiative losses (along with gravity 
and radiative driving terms), and directly relates to the work 
presented here. While those simulations have already shown the 
complex morphology of the circumstellar medium, we augment 
those results here for a hydrogen rich Wolf-Rayet (WNL)+0 bi- 
nary and a Luminous Blue Variable (LBV)+0 binary, paying 
specific attention to the role and nature of the occurring thin- 
shell instabilities. 

In this paper we focus on the circumstellar medium of mas- 
sive binaries with a relatively wide orbit taking the orbital pe- 
riod equal to 1 year. For the components of the binary system 
we choose a combination of a WNL and an O-type star and a 
combination of a LBV and an O-star, using generic parameters 
rather than attempting to simulate existing binary systems. To fa- 
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cilitate a comparison between the two systems, we use identical 
stellar mass and orbital parameters for both. This can be inter- 
preted as two stages in the evolution of a massive binary, with the 
most massive of the two stars making the transitio n from hydro- 
gen rich Wolf-Rayet to Luminous Blue Variable (lLanger et al.l 
Il994b . 

In the region where the winds collide, a shell forms, which 
is strongly compressed by the collision. For the WNL + O star 
binary the collision is adiabatic on both sides resulting in a thick, 
smooth shell. However, for the second binary, owing to the dis- 
parate nature of the winds (the LBV wind is relatively slow and 
quite dense, whereas the O-star wind is ten times faster, but has 
comparatively low density), the nature of the shock is radiative 
on one side and at least partially adiabatic on the other. Since 
the resulting shell is thin, as a result of efficient radiative cool- 
ing, it will be subject to combined hydro and radiative instabil- 
ities, depending on the (evolving) wind parameters. We identify 
the nature of the instabilities in the specific LBV+O scenario, 
which form in the shell due to both the wind interaction and the 
stellar orbital motion. In particular, the shell is characterized by 
the presence of multiple thin- shell instabilities, with clear differ- 
ences between leading and trailing interaction fronts. 

2. Generic binary models 

2.1. Governing equations 

We simulate the hydrodynamical interaction between the winds 
of two massive stars, solving for the dynamics in the center 
of mass frame of the two stars. The stellar masses, orbit and 
wind parameters are shown in table [T] The O-star and WNL 
win d mass-loss rate are chosen rather low , follow i ng the findings 
oflBouret etail (120051) : I Vink & de Koterl (l2005h : iMokiem et all 
d2007h and others that clumping has led to overestimation of the 
mass-loss rat e of hot stars. The LBV mass-loss rate reflects val- 
ues found bv lVink&deKoterl (120021) . 

Since we simulate a binary with a relatively wide orbit, we 
can ig nore the acceleration zone o f the stellar winds (typically < 
10R* lLamers & Cassine lli 1999, and citations therein) and 
model wind- wind interactions between winds at terminal veloc- 
ity. This also justifies the neglect of gravity and radiation -driven 
body forces, which were taken along in the work by Pittard 
(2009). The changing orbital positions are calculated following 
Kepler's law and the terminal wind zones are prescribed as ex- 
plained below. We take into account the effect of optically thin 
radia tive cooling, using the cooli ng curve for solar metallicity 
from iMellema & Lunqvisi (120021) . Since the medium close to 
two massive stars can be assumed to be photo-ionized, we keep 
the gas at a minimum temperature of 10000 K throughout the 
entire computational domain. 

We solve the usual equations of hydrodynamics: conserva- 
tion of mass, 

dp 

-£ + V • (pv) =0, (1) 
ot 

with p the mass density and v the velocity vector; conservation 
of momentum, 

^ + V.(pvv) = -Vp, (2) 
ot 

with p the thermal pressure, and with the energy equation given 
by 

d JL + v.(ev) + V-(pv) = -n 2 A(Tl (3) 



with e the total energy density, and energy loss due to radiation 
in the righthand side. This term depends quadratically on the ion 
density n, and A(T) is t he temperature-dependent c ooling rate 
for solar metallicity from|Mellema & Lunqvis t ( 2002), which is 
simila r to the SPEX-code based cooling curve used by [Pittard 
( 2009), but incorporates fewer ion species. 

2.2. Numerical setup 

We use the MPI-AMRVAC code dMeliani et al.l 120071) . using 
its hydrodynamical module, on grids where the local resolu- 
tion can be changed through adaptive mesh refinement (AMR). 
We implemented optically thin ra diative co o ling u sing the new 
exact integration method from iTownsendl (120091) . which im- 
proves calculation speed and numer ical stability. A compari- 
son by Ivan Marie & Keppensl (120101) between different source 
term treatments for these optically thin radiative loss terms, in 
ID to 2D single wind bubble evolutions, clearly showed the 
need for AMR combined with robust radiative cooling prescrip- 
tions. For the hydro, we use a shock-capturing TVDLF scheme 
dToth & Odstrcill Il996l) . employing limited linear reconstruc- 
tions. 

Our grid consists of a flat slab measuring 2.5 x 10 14 cm along 
the X- and Y-axes (the orbital plane of the binary) with 240 
grid points along each axis in the crudest refinement level and 
2.5 x 10 13 cm along the Z-axis with 20 grid points at its low- 
est resolution. For this particular simulation this means that in 
the X-Y plane the grid covers approximately 4 times the orbital 
separation along each axis. Assuming a typical O-star radius of 
about 10R o , this amounts to 350Ro- s tar- The radius of WNL 
and LBV stars is less well-known, as these stars have optically 
thick winds. Therefore, the effective temperature is not defined at 
the stellar surface, but somewhere in the wind region. Typically 
though, the WNL radius will be somewhat larger than the ra- 
dius of an O-star while the LBV radius will be about an order of 
magnitude larger than Ro-star- 

We then use two additional grid levels, where the intermedi- 
ate level achieves an effective resolution of 480x480x40, which 
is controlled by the automated error estimation procedure based 
on density variation, while the top grid level is at all times en- 
forced in the direct vicinity of the two stars, making the effective 
resolution there 960 x 960 x 80. The latter grid refinement level is 
thus activated in a fully user-controlled manner, employing the 
knowledge of the instantaneous stellar positions. 

At the beginning of each time step, this orbital position is 
calculated from the orbital parameters, and a sphere around each 
star of 5 x 10 12 cm in diameter is filled with free- streamin g 
wind material. This approach mimics the work by Pittard ( 2009), 
while external to this radius, the wind expands freely according 
to the governing hydrodynamics. 

3. Circumbinary medium morphology 

3.1. WNL+O binary 

The morphology of the circumbinary medium of the WNL+O 
binary, as shown in Figs.Q]and[2]is a typical exampl e of a purely 
adiabatic wind- wind collision (IStevens et al.lll992h . The shock 
is adiabatic on both sides of the contact discontinuity, resulting 
in a thick shell comprised of shocked wind gas. The edges of the 
shell are defined by the point where the thermal pressure inside 
the shocked wind region is equal to the ram pressure of the wind. 
The WNL wind, which has a much higher ram pressure (a factor 
of 7.5) clearly dominates the collision and has folded the O-star 
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Table 1. Binary parameters. All parameters regarding the primary star denoted with subscript 1, all parameters regarding secondary 
with subscript 2. 



Simulation 


Massi 


Mass 2 


Mi 


M 2 




V2 


Period 


eccentricity 




[Mo] 


[Mo] 


[Me/yr] 


[M /;yr] 


[km/s] 


[km/s] 


[yr] 




WNL+O 


50 


30 


5 X 10" 6 


5 x 10" 7 


1500 


2000 


1 





LBV+O 


50 


30 


1 x 10" 4 


5 x 10" 7 


200 


2000 


1 







Fig. 1. Gas density in g/cm 3 of the WNL+O binary system af- 
ter one full orbital revolution (lyr). Owing to the higher ram 
pressure of the WNL (violet sphere) wind the bowshock curves 
around the O-star (white sphere). The shell is thick on both sides 
of the contact discontinuity and shows no sign of instabilities. 




Fig. 2. Thermal pressure in erg/cm 3 for the WNL+O binary sys- 
tem at the same time as Fig. [T] The thermal pressure of the 
shocked wind region balances against the ram pressure of the 
two winds. 



wind back into a bow-shock. Owing to its thickness, the shell is 
not subject to thin- shell instabilities and remains smooth. 

The leading edge of the shell (lower part in Figs. [T] and [2]) 
runs approximately parallel to the WNL wind, whereas the trail- 
ing edge falls slightly behind due to the orbital motion, which 
runs counter-clockwise. Since there is a shear along the edges 
of the shell it would, in theory, be subject to Kelvin-Helmholtz 
instabilities. However, the probability of such instabilities oc- 
curring is reduced by the orbital motion of the shell, which is 
perpendicular to its surface and tends to wipe out such instabili- 
ties. 



Fig. 3. 3D image of the circumstellar medium around an LBV 
+ O binary, after one orbital period (1 yr). The slice along the 
equatorial plane shows the logarithm of the density. An isosur- 
face along which the absolute velocity equals 300 km/s is repre- 
sentative of conditions along the contact discontinuity between 
the shocked LBV wind and the shocked O-star wind. N.B. the 
z-axis has been stretched by a factor 2 to show more clearly the 
3D nature of the instabilities. 



3.2. LBV+O binary 

The morphology of the circumstellar medium around the 
LBV+O binary system is presented in Figs. [3] through [6] and 
clearly shows a much more complicated structure due to insta- 
bilities in the shell. Figure[3]gives an impression of the 3D struc- 
ture through a constant absolute velocity isosurface, at a value 
that selects the contact discontinuity between the two shocked 
winds. Since this follows the contour of the shell, it shows the in- 
trinsic three-dimensional structure of the instabilities. Figures [H 
through [6] show 2D slices in the orbital plane, taken after a full 
orbital evolution. The LBV wind, which has the higher ram- 
pressure ( J P ra m(LBV)/ J P ra m(0) = 20), dominates the interaction, 
so a bow- shock forms ahead of the O-star as it plows into the 
LBV wind, but unlike the WNL+O binary where this bow-shock 
looks almost symmetric, the LBV+O binary shows that the trail- 
ing edge of the bow- shock has fallen behind in the orbital mo- 
tion, despite the fact that the LBV wind has a higher ram pressure 
than the WNL star from the first simulation. This is the result of 
the low velocity of the LBV wind and will be discussed in more 
detail in Sectionl4~2l 

Figure |7] demonstrates the evolution of the circumstellar 
medium over time as the LBV+O binary completes its sec- 
ond orbit, showing the density of the gas after 1.5 (left) and 2 
(right) orbits respectively. The general morphology of the shell 
remains the same, though the thin-shell instabilities at the front 
of the bow-shock become a bit more pronounced. This is fur- 
ther demonstrated in Fig.[8l which shows the development of the 
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Fig. 4. Density of the circumstellar medium in the LBV + O sys- 
tem at the same moment in time as Fig. [3] This figure shows a 
cross-section in the orbital plane of the binary. The wind from 
the Luminous Blue Variable (blue sphere) dominates due to its 
higher ram-pressure. The shell created by the collision is unsta- 
ble, showing two kinds of thin shell instabilities: non-linear or- 
der thin-shell, instabilities at the head of the bow shock around 
the O-star (white sphere) and in the region directly between the 
two stars. Downstream (red arrow) the instabilities are of lin- 
ear thin- shell type. Even further downstream, after the shell has 
made its curve around the approaching LBV star the shell begins 
to diffuse due to its own internal pressure. 

LBV+O binary shell over time. It presents isosurfaces defined 
at an absolute velocity of 300km/s (see also Fig. [3j at 5 sub- 
sequent moments in time, covering one complete orbit. The first 
contour (red) shows the shell 1/5 of an orbit after Figs. [3] through 
[6|and the final contour (white) shows the shell at the same time as 
Fig. [TJ The global shape of the shell remains constant over time, 
rotates rigidly around the center of mass of the binary. The local 
instabilities vary over time, but maintain the same general shape, 
with relatively small deformations in the region directly between 
the two stars and much larger in the leading and trailing parts. At 
the leading end of the shell the deformations also maintain their 
'saw-tooth' appearance. The distortions of the shell directly be- 
tween the two stars remain comparatively small. This is due to 
two effects: As the distortion moves toward either of the stars it 
will experience an increasing ram pressure owing to the increase 
in wind density, which slows down its motion. Also, the growth 
time of these instabilities is on the order of the orbital period, so 
the shell itself moves significantly before the instability has time 
to grow further. The nature of these instabilities is explored in 
more detail in Section 

4. The collision region 

4.1. Nature of the shock 

As shown by lStevens et al.l (Q992), the structure of the interac- 
tion region is determined by the nature of the shocks. For an adi- 
abatic shock an extended shocked wind region with a high tem- 
perature will develop between the free- streaming wind and the 
contact discontinuity. On the other hand, if the shock is radiative, 
the shocked wind region will be compressed into a thin shell. 
The defining crite rion for colliding wind shocks is the cooling 
parameter, which IStevens et al.l ([l992) determined as the rela- 
tion between the cooling timescale t coo \ = kT s /(n w A(T s )) and 




Fig. 5. Temperature of the circumstellar medium in K at the time 
from Figs. [3] and IU This figure shows the reason for the split in 
thin- shell instabilities. At the front of the bowshock as well as 
between the two stars, the thermalized zone is extremely thin, 
as a result of effective radiative cooling. The shock, owing to 
its high density and temperature is almost completely radiative 
there, making the shell subject to ram-pressure on both sides. In 
the trailing end (red arrow) the shock interaction is less violent 
and densities are lower. As a result the shock is no longer ra- 
diative and a relatively thick thermalized zone forms. Hence the 
zone between contact and termination shock for the LBV wind 
bubble feels ram-pressure on one side and thermal pressure on 
the other. This changes the nature of the thin- shell instabilities. 




Fig. 6. Thermal pressure in erg/cm 3 of the circumstellar medium 
at the same time as Figs. [3] through This plot shows how the 
thermal pressure supporting the shell from the O-star side dis- 
appears at a larger distance (blue arrow). The shell still feels the 
ram pressure from the LBV star wind, which pushes it ahead, but 
without a counter force from the other side it is no longer com- 
pressed and starts to expand under its own internal pressure. This 
plot also shows how the O-star wind forms bowshocks around 
the individual deformations in the downstream region. 



the escape time from the shocked region t QSC = d/c s , with k 
the Boltzmann constant, T s the shocked wind temperature, n w 
the particle density in the wind at the shock, d the distance 
to the star and c s the post-shock sound speed. Assuming that 
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c s ~ v ~ ^/T~ s and taking A constant for typical shock tempera- 
tures of 10 7 ...10 8 K, this gives us a cooling parameter of: 



with vs the wind velocity in 10 8 cm/s, dn the distance between 
star and shock in 10 12 cm and M_ 7 the wind mass-loss rate in 
10~ 7 M o /yr. Although this equation is not very precise (it ne- 
glects the temperature dependence of the radiative cooling), it 
gives an indication as to the nature of the shock. 

For the WNL+O binary, the cooling parameter for the O- 
star is: x ~ 50 in the region directly between the two stars, 
whereas for the WNL star it v&x ~ 4, indicating that both shocks 
are at least partially adiabatic, as x > 1- This is confirmed by 
the compression factor, which for both shocks is approximately 
4. Obviously, further away from the collision point the shocks 
becomes even more adiabatic as the distance to the star increases 
leading to lower densities and less efficient cooling. 

For the second binary, directly between the two stars the 
LBV wind has x ~ 10~ 4 ^ 1> which makes this shock com- 
pletely radiative, resulting in a thin shell of stalled LBV wind 
material. The O-star wind again has^ ~ 50 > 1, which makes 
the shock adiabatic. Compression rates vary due to the local in- 
stabilities. However, though the O-star wind does develop a ther- 
malized layer between the free- streaming wind and the contact 
discontinuity, this layer is comparatively thin (approx. twice the 
cross-section of the compressed LBV wind). This is of no great 
consequence for the WNL+O binary, which doesn't show insta- 
bilities, but is crucial in the LBV+O binary. 

This calculation does not take into account the orbital mo- 
tion of the stars, which greatly complicates the situation for the 
LBV+O binary, where the slower of the two winds (the LBV 
wind) moves with a speed that is comparable to the orbital ve- 
locity of the rigidly rotating bowshock. For the WNL+O binary 
the shell is effectively caught between the ram pressure of both 
winds along its entire length. However, for the LBV +0 binary 
the situation is completely different. At the edge of the bow- 
shock, where it hits the LBV wind material (white arrow in 
Figs. IH through O, the shell is bent further toward the O-star 
due to the higher ram pressure of the LBV wind as well as the 
fact that it is sweeping up material that is much denser than in 
the case of the WNL wind. As a result, the relevant collision ve- 
locity on the LBV wind side of the contact discontinuity is not 




Fig. 8. The shell between the LBV and O star at five different 
points in time, starting l/5 th orbit after Figs. [3] through (red) 
and at each l/5 th orbit afterwards (green, blue, yellow) The final 
shell (white) coincides with the righthand plot in Fig. [7J The 
shells are defined as an isosurface contour drawn at a velocity of 
300km/s (see also Fig. [3]). The global shape of the shell remains 
constant over time, rotating rigidly around the center of mass 
of the binary (0,0). Local instabilities vary, but retain the same 
general characteristics, being small between the stars, but much 
larger in the advancing and trailing ends of the shell. 

the LBV wind speed, which runs almost parallel to the shock 
front, but the orbital velocity of the shock into the wind mate- 
rial. At the trailing end of the bow shock (red and blue arrows 
in Figs. IH through [6]) , the shell has to keep up with the orbital 
rotation while it is being driven ahead by the slower of the two 
winds (either the WNL or LBV wind). This can result in the loss 
of the shock as will be discussed in Section l4~2l 

4.2. The loss of ram pressure balance in the trailing part of 
the shell 

The shape of the trailing end of the collision region depends on 
the relative velocities of orbital motion and the stronger of the 
two winds. This wind drives the trailing end of the bow shock 
ahead, which means that the shell can never move faster than the 
component of the LBV wind velocity perpendicular to the shell 
(or even l/4 th of this velocity as long as a shock is maintained 
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at the collision front, because of the Rankine-Hugoniot condi- 
tions). Since the shell is, effectively, rigidly rotating around the 
center of mass of the binary (see also Fig. [5]), its orbital velocity 
has to increase linearly with the distance to the center of mass 
(v orbit - 2nRQ). Therefore, for larger radius the orbital veloc- 
ity approaches the LBV wind velocity. At this point the shell 
has to fall back to catch a larger percentage of the LBV wind 
momentum (which would otherwise hit it at a very shallow an- 
gle). Eventually even this fails and the shell starts to lag behind 
the orbital motion. A similar eff ect can be observed in some of 
the simulations bv lPittardl (l2009h . The WNL wind, which is 7.5 
times faster than the LBV wind, only reaches this point at the 
edge of our computational domain, which can be seen in the 
gradual change in the curvature of the shell in Figs. Q] and [2] at 
a distance of about 10 14 cm from the center of ma ss. This agrees 
with the description by Parkin & Pittardl (120081) . which shows 
that the angle of the bow- shock with the axis between the two 
stars depends on the velocity of the slowest of the stellar winds, 
relative to the orbital motion. 

The WNL+O binary maintains ram pressure balance be- 
tween the two winds in the entire computational domain. 
For the LBV+O binary the situation is completely different. 
Downstream from the collision (red arrow in Figs. |H through [6]), 
the shocked wind region effectively no longer feels the ram pres- 
sure of the O-star wind (which runs parallel to the shock), allow- 
ing the shocked wind region to expand. We can set limits for the 
shell displacement. If we assume that enough thermal pressure 
always remains to counterbalance the ram pressure of the LBV 
wind, then the shell stays at a constant distance to the LBV star 
and its displacement equals the motion of the LBV star perpen- 
dicular to the shell. This is the absolute minimum displacement. 
The maximum can be deduced from a hypothetical case in which 
there is no thermal pressure at all and the shell is simply carried 
along by the LBV wind. In that case the displacement is equal 
to the LBV wind velocity times the time interval. Obviously this 
overestimates the motion, since the mass of the shell itself (and 
its related inertia) would slow it down even without any thermal 
pressure. Using these two estimates we obtain for a quarter orbit 
(90°, see Fig. 0), 

Atf min = tf orblt (LBV) = 1.75 x 10 13 cm, (5) 
Atf max = i;(LB V) A t = 1.58 x 10 14 cm. (6) 

Clearly, these two limits lie far apart and the actual displacement 
of the shell in our simulation, which is about 5 x 10 13 cm relative 
to the center of mass of the binary (point (0,0) in Fig. [9]), falls in 
between. The Rankine-Hugoniot conditions determine that the 
post shock speed is l/4 th of the pre-shock speed. Therefore, the 
shell could only move l/4 th of AR max if a shock was maintained 
between the shell and the LBV wind. This would limit the dis- 
placement of the shell to AR = 3.95 x 10 13 cm. Since the actual 
displacement is larger, we can conclude that the shock between 
the LBV wind and the shell has been lost. 

The point where the shell will deviate from the balance be- 
tween the ram pressure of the two winds can be approximated by 
using the condition that the component of the LBV wind along 
the orbital direction of motion has to exceed the orbital velocity 
determined by the binary stars: v L > 2nR£l, with v L the com- 
ponent of the LBV wind along the direction of orbital motion. 
This condition has been worked out in Appendix |A) For our 
LBV+O binary the resulting equation is easy to use, since the 
shell between the two stars runs almost perpendicular to the axis 
connecting the two stars. Using an approximation where AX in 
Eq. IA.5l is constant (set to 2 x 10 13 , the distance from the cen- 
ter of mass to the point where the two winds collide head-on), 



the turning point of the shell lies at A F ^ 2.35x 10 13 cm. This 
coincides quite well with the point where our simulation shows 
a significant deviation from the bowshock in Figs. |4] through [7] 
From this point on the shell starts to curve backwards into the 
LBV wind region until, eventually, it runs parallel to the orbital 
motion. 

This calculation does not take into account the distortion of 
the shell due to instabilities. Even a comparatively small insta- 
bility can effectively shield a large section of the shell from the 
ram pressure of the O-star wind as a bow-shock forms around 
the instability as can be seen in Figs. |4] through H 

As the shell curves away from the O-star wind (blue arrow 
in Figs. |4] through [6]) , because it cannot keep up with the orbital 
rotation, it loses the support from the thermalized O-star wind 
region and starts to expand under its own internal pressure. Since 
the O-star wind moves faster than the shell, the area behind the 
shell is essentially empty making it easy for the shell to expand. 
It gets pushed ahead by the wind of the approaching LBV star 
and accelerates away as the ram pressure of the wind overcomes 
the inertia of the shell. For the WNL+O binary, this would only 
happen at a much larger distance from the binary stars, as the 
faster WNL wind is capable of maintaining a shock at higher 
orbital velocity. 

4.3. Ionization state 

Our assumption of complete photo-ionization is supported by 
the results. The particle density in the shell reaches a maxi- 
mum of about 10 10 cm~ 3 . For a typical O-star, which produces 
on the order of 10 49 s" 1 high energy photons, this would im- 
ply a Stromgren radius of about 5 x 10 13 cm (Eq. 5.14 from 
Dys on & Williamslll997l) if the density was constant from the 
O-star to the outer edge of the shell. In reality the density around 
the O-star is much lower due to the free expansion region of 
the wind, leaving only the thin shell to absorb the photons 
sinc e recombination in the free- stre aming wind region is negligi- 
ble dGarcia-Seg ura & Francoll 1996b . This means that the O-star 
alone can ionize the shell between the two stars, even without the 
presence of the LBV star. Only at larger distances from the star 
(once the matter has been carried away downstream) will it re- 
combine. Obviously, this is a generalization. Inside the shell high 
density clumps may form, which will be able to recombine due 
to self- shielding if they become optically thick. However, such 
structures will be small and can not be resolved by our current 
simulation. Moreover, this calculation fails to take into account 
the presence of dust which can shield the gas from the ionizing 
radiation. 



5. Multiple instability characterization for the 
LBV+O binary 

For the WNL+O binary both shocks are (near-)adiabatic, lead- 
ing to a relatively thick shell that shows no instabilities. The 
LBV+O binary on the other hand reveals a complicated struc- 
ture, as the highly compressed LBV wind forms a thin she ll that 
is subject to thin-shell instabilities (iVishniacll 19831 [19941) . That 
shells between co lliding stellar winds a re prone to these instabil- 
ities was found by lStevens et al.l {1992) in pioneering 2D simula- 
tions, which did n ot fully account for the effect of orbital motion. 
3D simulations by lDgani et al.l (Il996l) showed similar results for 
the collision between a stellar wind and a parallel flow, as can 
occur if the star itself moves with supersonic speed relative to 
the interstellar medium. 
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Fig. 9. Displacement of the shell over a 1/4 1 orbit (0.25 yr). The 
distance A R over which the shell moves is determined by its own 
inertia, the remnant of thermal pressure that balances it against 
the ram pressure of the approaching LBV star and the velocity 
of the LBV wind, placing an upper limit on the shell velocity. 



Our 3D simulation of two stars in orbit shows a more com- 
plicated pict ure, also foun d in some of the double O-star models 
analyzed by lPittardl (120091) . For our LBV + O star system in cir- 
cular orbit the shell ins tabilities come in two forms: non-linear 
thin-shell instabilities (Vishniac 1994), which occur if a thin 
shell is caught betwee n ram-pressure f rom both sides, and linear 
thin- shell instabilities (iVishniacll 19831) . which depend on having 
ram-pressure on one side and thermal pressure on the other. As 
the winds from both stars are relatively cool (~ 10 000 K), the 
only way to create hot gas is if it gets thermalized by a shock. 
Whether such a thermalized zone forms, depends on the effec- 
tiveness of the radiative cooling, as discussed in Sect. [3] 

Directly between the two stars the winds collide head on. 
The shock is completely radiative on the LBV side and at most 
marginally adiabatic on the O-star side. As a result the shocked 
LBV wind is compressed to a very thin shell (estimated com- 
pression factor « 10. ..20). The shocked O-star wind is also 
compressed, but not to the same extent (estimated compression 
factor « 5... 10). Initially, the instability is of the linear thin- 
shell variety as the shocked LBV wind shell is caught between 
the LBV wind ram-pressure and the thermal pressure from the 
shocked O-star wind. However, becau se of the long orbital pe- 
riod the instabilities have time to grow (Vishniac 1983, eqn. 2.23 
gives a growth rate of about 1 month for this specific case, com- 
pared to the 1 yr orbital period.). As the instabilities cross the 
thin shocked O-star wind layer (~ 5x 10 12 cm), the reverse shock 
in the O-star wind will start to conform to their shape. At this 
moment the nature of the instabilities becomes non-linear as they 
now effectively feel ram-pressure from both sides. The growth 
rate of such instabilities is no longer determined by the sound 
speed in the shell, but by the velocity of the wind (200 km/s for 
the LBV wind) which is faster than both the sound speed and the 
orbital motion of this part of the shell (< 50 km/s), allowing the 
instabilities to keep growing. 

At the front of the bow shock (white arrow in Figs.[4]through 
[6]) the shell does not actually feel the LBV wind, which moves 
nearly parallel to its surface. Instead, it feels a ram-pressure 
caused by its own movement (at an orbital velocity of about 
90 km/s) into the LBV wind region. However, as the develop- 
ment of the instability causes deformations that extend into the 
free- streaming LBV wind region, they are subjected to a side- 
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ways force by the ram-pressure of the LBV wind, causing a 
strongly asymmetric shape, which reflects the relative forces act- 
ing on them. In the largest deformation, near the white arrow in 
Figs. |H and [6l the slope facing the LBV wind is approximately 
four times as long as its displacement from the shell, conforming 
to the 4:1 ratio between the LBV wind ram-pressure and the or- 
bital motion ram-pressure (vLBVwind ~ 2y or bit an d Pram ~ v 2 .) The 
other instabilities in this region are less asymmetric as the first 
instability shields them at least partially from the LBV wind. 
These instabilities look somewhat similar to Kelvin-Helmholtz 
instabilities, but a close-up of this same region at the same time 
as Fig. [7] (right) shows that there is no circular motion in them 
(Fig. [TOb . This figure shows streamtraces of the velocity in the 
equatorial plane. The wind velocity clearly dominates but is de- 
flected by the local instabilities. The interface between the shell 
and the LBV wind meets the criteria for Kelvin-Helmholtz in- 
stabilities due to the shear between the shell and the wind, but 
the shell moves perpendicular to the shear interface, effectively 
crushing such instabilities even if they start to form. Owing to the 
bow-shocks that form around the instabilities on the O-star side 
of the shell, pockets of hot gas form behind the individual de- 
formations, which therefore start to feel thermal pressure, rather 
than purely ram pressure. This adds a linear thin-shell instability 
component. As a result the typical wavelengths here are longer 
than in th e purely non-linear thin shell instabilities between the 
two stars (IVishniaclll9"83ill994l:[Garcia-Segura et alJl!996bb . 

At the trailing end of the shock (red arrow in Figs. |H through 
[6]) the thermalized O-star wind forms an extensive layer, so that 
the shell only feels ram-pressure from the LBV wind, counter- 
balanced by thermal pressure from the shocked O-star wind. As 
a result, the instabilities here are becoming linear thin- shell in- 
stabilities, leading to longer wavelengths of the individual de- 
formations of the shell. As these instabilities extend into the 
shocked O-star wind, they feel a shear force due to the velocity 
within the shocked O-star wind region, which runs parallel to the 
shell. This causes the beginning of Kelvin-Helmholtz type insta- 
bilities. In this region the temperature of the shell could drop be- 
low our own limit of 10000 K as the shocks here are relatively 
weak and local clumps may become dense enough for recom- 
bination. In that case it would be a good environment for dust 
formation as the local density is still high, but the temperature 
could drop to the order of a f ew thousand K, which would be 
acceptable for dust formation (ICherchneff et al.|[2000b . Around 
the individual instabilities the O-star wind forms local bow- 
shocks. This effectively decreases the thermal pressure on the 
shell, which starts to expand to the point where it can no longer 
be considered thin. It might be possible for individual clumps of 
dense matter from the shell to break off due to the shear force 
from the O-star wind and disperse into the low density region. 
However, this is difficult to predict from numerical simulations, 
since the resolution of the grid effectively determines the mini- 
mum size of a high density element. 

Even further downwind (blue arrow in Figs. [4l through© the 
shell curves away from the O-star wind due to the orbital motion 
of the stars. At this point it no longer feels any sort of pressure 
from the O-star wind and therefore is no longer subject to first 
order thin- shell instabilities. The region behind the shell (where 
the O-star wind used to be) has a very low density, since the O- 
star wind is much faster than the motion of the shell and leaves 
a rarefied area into which the shell moves. It is being pushed 
ahead by the wind from the approaching LBV star, but with- 
out supporting pressure on the other side it is no longer being 
compressed and starts to diffuse as its own internal thermal pres- 
sure (primarily caused by the high density) causes it to expand. 
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Fig. 10. Streamtraces of the gas velocity in the orbital plane at 
the same time as in the right panel of Fig. H The instabilities 
deflect the LBV wind. There is no sign of Kelvin-Helmholtz in- 
stabilities (no circular motion). Instead, the asymmetric shape of 
the instabilities is caused by the shear of the LBV wind running 
parallel to the shell. 



Since the high density parts have more inertia than the low den- 
sity regions, they will tend to lag behind the low density regions, 
causing further distortion of the shell. At this point the transition 
between the LBV wind and the shell is no longer a shock, so the 
shell has effectively become a high density ridge, moving ahead 
of the wind. 



6. Conclusion and outlook 

Our simulation shows that the shell created by the collision of 
two different stellar winds can be subject to at least two kinds 
of thin- shell instabilities simultaneously, depending on the wind 
and orbital motion parameters. We find that the orbital motion 
plays a crucial part in shaping the interaction region between 
two stars, both on a global scale, in determining the position of 
the shell, and on a local scale, in determining the nature of the 
instabilities. 

6.1. Observational properties 

In terms of observations, the LBV+O binary would primarily 
radiate in the optical/IR. The (at least partially) adiabatic nature 
of the fast wind shock, will limit the production of X-rays. The 
radiatively cooling LBV wind shock would show a strong signa- 
ture in the visual part of the spectrum (and possibly UV, though 
again, the high density of the LBV wind could absorb most of 
this). Due to the stability of the general shape of the interaction 
region, the optical signal, as produced by the shock, would prob- 
ably be quite stable. Any strong variations would be more likely 
to be due to absorption by optically thick parts of the shell or 
LBV wind as they pass between the shock and the observer. The 
dense, rapidly cooling shell can also enable dust formation, as 
observed in many massive binary systems. This le ads to a clear 
IR signature as can be obs erved in WR 140 (e.g. lMonnier et al.l 
2002; V arricatt et al.ll2004l) and the far more massive Eta Carinae 
dSmithl 12010). The emission in optical will be increased due to 
the presence of extensive instabilities, which cause more kinetic 
energy to be transferred to thermal energy, which in turn radi- 
ates away through cooling. Since the instabilities enhance the 



local density in the shell, they will also increase the formation of 
dust, which in turn increases the IR radiation. 

The WNL+O binary may produce more high energy pho- 
tons as the fast WNL wind shock, though technically adiabatic, 
comes close to being radiative in the region directly between the 
two stars. Whether any such photons would be observable is an- 
other matter. They have to penetrate either the high density shell, 
or the free- streaming region of the WNL wind, which is at least 
partially optically thick. The result, if any X-ray emission would 
be seen at all, would depend very much on the angle of the ob- 
server with the shock cone and could lead to a periodic signal 
such as observed for Eta Carinae (ICorcoranl2005l: lOkazaki et all 
120081) . though in the case of Eta Carinae r adiative braking o f the 
wind may play an important role as well dParkin et al .1120091) in 
determining the number of photons produced. This is not an is- 
sue in our binary where the winds only interact at terminal veloc- 
ity. The generally adiabatic nature of the shocks in the WNL+O 
binary prohibits the production of a strong optical signal and 
the luminosity of the massive stars themselves will tend to over- 
whelm any sign of the shock in the optical, unless the binary 
system could be spatially resolved. The most promising way to 
observe the shock would again be in the IR band, tracking the 
dust formed in the shell, though the amount of dust would be 
much less than for the LBV+O binary due to the lower densities 
in the shell and the fact that the WNL wind would contain less 
carbon than the LBV wind. 

6.2. Future work 

We here expand on the 2-D findings by IStevens et al.l (Il992l) . 
showing thin- shell instabilities in 3-D simulations. Obviously, 
the exact nature of instabilities will depend on the character of 
the colliding winds and can be expected to change if different ini- 
tial parameters are used. However, our results make clear that the 
orbital dynamics of the binary, which cause the split in the na- 
ture of the thin-shell instabilities for the LBV+O system, cannot 
be ignored. Future work will consist of performing more para- 
metric studies of the occurrence and dynamical influence of such 
multiple instabilities. 
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A.1. Mathematical analysis 

The exact point where the trailing end of the shell between the 
two stars starts to bend away from the ram pressure balance is 
a function of the binary parameters. Its location depends on the 
velocity of the stronger of the two winds, which is pushing it 
ahead in the orbit, and how it scales against the orbital velocity. 

We will use a 2-D coordinate system with the two stars po- 
sitioned on the X-axis, rotating counter clockwise, and the cen- 
ter of mass of the binary system in the origin (see Fig. lA.lb . 
Close to the binary stars, the location of the shell can be de- 
scribed as a collection of points (AX, A Y), relative to the cen- 
ter of mass of the binary, which follow the curve along which 
the ram pressure of the two winds is equal. However, as was 
already mentioned in Sect. [3] there is an upper limit to the or- 
bital velocity of the shell: it cannot exceed the velocity of the 
wind that pushes it. (In reality it will always move slower due 
to the inertia of the shell.) Consequently, the shell will start to 
deviate from the ram pressure balance curve at the point where: 
Vi < 2tiR£1, with v L the component of the wind velocity along 
the direction of orbital motion, Q the angular velocity of the bi- 
nary and R = V A X 2 + A Y 2 the distance from the center of mass 
to that particular point on the shell. 

This is shown in Fig. IA.U Here the orbital motion in point 
(A X, A Y) and the wind velocity (along the line from the LBV 
star to the point (AX, A Y)), make a steep angle, indicating that 
only a small part of the wind momentum actually pushes the 
shell ahead in orbit. The exact fraction follows from the binary 
parameters. If a is the angle of the wind with the horizontal 
line between the two stars and J3 is the angle of the position of 
(AX, A Y) with the center of mass (0, 0), the angle of the wind 
with the orbital motion is given by 
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where R(LBV) = V A Y 2 + (AX + X(LBV)) 2 is the distance 
from the LBV star to the point on the shell. The component of 
the LBV wind along the direction of orbital motion, v L is now 
given by: 
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Appendix A: Turning point of the shell 

For binaries where one of the winds is (relatively) slow com- 
pared to the orbital motion, the location of the interaction front 
between the two winds can be determined more by the relative 
velocity of this slow wind, than by the point where ram pres- 
sure balance exists between the two winds. Here we explore the 
situation where the slower wind is also the stronger in terms of 
mechanical luminosity. 



Since R = VAX 2 + A Y 2 , this equation can be evaluated for any 
given position along the X-axis to yield the value A Y at which 
the wind can no longer support the orbital motion. 

The important parameters here are the relative velocities of 
the rotation and the wind as well as the relative masses of the 
stars. In the most extreme case, where the LBV mass completely 
dominates the binary, R(LBV) = R so the right hand term of 
Eq. IA.5I goes to zero and the wind would never be able to push 
the shell ahead in orbit. If the LBV actually has a smaller mass 
than its companion, R(LBV) » R, the region where the wind 
can push the shell into orbit becomes large. 
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Fig.A.2. Rotational velocity (continuous) and LBV wind 
(dashed) along the orbital direction of motion for AX = 2 x 
10 13 cm in the LBV+O binary. The curves cross each other 
twice, showing the roots of Eq. IA.5I The second root (AY ^ 
2.35 x 10 13 cm) marks the point where the collision shell starts to 
deviate from a ram pressure induced bowshock. 



A.2. Application to a binary simulation 

The condition in Eq. IA.5I can be used to determine the point 
where the shell between the components of a binary system will 
deviate from the ram pressure balance between the two stellar 
winds. 

Let us assume, for the moment, that the position of the stars 
is stationary. This is a reasonable assumption for a binary system 
where the wind velocities are much larger than the orbital ve- 
locity. Under these conditions the location of the collision front 
between the stellar winds can be found analytically as a collec- 
tion of points (A X, A Y) where the ram pressure of the two winds 
is equal. This curve will follow a bow shock that is symmetric 
around the X-axis. 

We can then use the coordinates of each point on this curve 
to find the point wher e this symmetrical bowshock no longer 
satisfies condition |A.5l At this point the shell will have to deviate 
significantly from the ram pressure balance, since it has reached 
its maximum velocity and can not keep up with the rigid rotation. 

As an example we use our simulation of an LBV+O binary. 
As shown in Figs. [4] through [7] the trailing end of the shell runs 
almost vertical between the two stars. Taking for example the X- 
axis coordinate of the collision point of the two winds directly 
between the two stars (AX ^ 2 x 10 13 cm) and the binary pa- 
rameters from Table Q] we find that the rotational velocity (left 
side of Eq. IA.5I) and the r otatio nal component of the LBV wind 
velocity (right side of Eq. IA.5I) follow the two curves shown in 
Fig. lA.2[ Since the equation is second order, there are two roots, 
one at about 1.20 x 10 13 cm and the second at 2.35 x 10 13 cm. 

Therefore, if we would follow the curve of ram pressure bal- 
ance we can determine three separate regions. Close to the axis 
con necti ng the two stars, A Y^eii is lower than the first root of 
Eq. lA.5l for the local value of AX. Here the wind is almost per- 
pendicular to the shell, whereas the orbital motion is almost par- 
allel to the shell. Obviously the wind cannot impart orbital mo- 
tion to the shell, which is effectively standing still, caught be- 
tween the two winds. Since both winds have a momentum com- 



ponent parallel to the shell and counter to the orbital motion (for 
A Y > 0) the matter is 'squeezed out' from between the stars 
until it reaches the point where (A X, A Y) S heii falls between the 
two roots of the equation. In this region the shell is caught up 
by the wind, which can push it into orbit. As the matter travels 
further downstream it reaches the second root, where the LBV 
wind becomes insufficient and it starts to fall behind the orbital 
movement of the binary. 

This analysis still ignores part of the orbital motion in that 
it assumes that the wind impacting on the shell at a given time 
t was also launched from the star at that same time t. In reality 
the wind was emitted earlier, which means that the star was at 
a different position. How important this is depends on the total 
displacement of the star between the time the wind is launched 
and the time it hits the shell. 



